F= 135; % Force in newtons. kg m/s^2
r = .036;   % Moment arm at which muscle is acting in distance from joint in meters.
theta = 0.0; % hanging down at side.
theta_offset = 5.0*pi/180; % muscle insertion point is not colinear with bone. Angle defined by bone and tendon.
g=9.8; % gravity in Newtons.
forearm_length = .33; % length in meters.
forearm_mass = 3; % assume it weighs 3kg.

theta_max = 130*pi/180;

tau = 1/100;
sim_seconds = 10;
theta_vel = 0;

theta1=zeros( 1, ceil(sim_seconds/tau) );

ts = tau:tau:sim_seconds;
figure(1);
for t = 1:ceil(sim_seconds/tau);

   tq_muscle = F*r*sin(theta+theta_offset);

   tq_load = g*sin(theta)*forearm_length*.5*forearm_mass

   tq_net = tq_muscle - tq_load;

   theta_accel = tq_net/(forearm_mass*(forearm_length*.5)^2);

   theta_vel = theta_vel + tau* theta_accel;
   theta = theta + tau* theta_vel;

   % limits
   if theta < 0
       theta = 0;
   end
   if theta > theta_max
       theta = theta_max;
   end
   
   [x y] = pol2cart(theta-pi/2, forearm_length);
   plot([0 x],[0 y])
   xlim([-forearm_length forearm_length]);
   ylim([-forearm_length forearm_length]);
   drawnow;
   pause(.01)

   % save
   theta1(t) = theta;

end

figure(2);
plot(ts,theta1*180/pi);